# =============================================================== #
# Reproducibility of:                                             #
# Drivers and Trends for the Equality of Opportunity for Sexual   #
# and Gender Minorities: A Panel Approach Equality of Opportunity #
# for Sexual and Gender Minorities 2024                           #
#                                                                 #
# Code written by Omar Alburqueque and reviewed by Paola Ballon   #
# Contact: oalburquequechav@worldbank.org, pballon@worldbank.org  #
# =============================================================== #

# ------- #
# Table 8 #
# ------- #

library(readxl)
library(tidyverse)
library(cluster)
library(e1071)
library(clusterCrit)
library(writexl)

project_root    <- getwd()
intermediate_dir <- file.path(project_root, "intermediate files")
outputs_dir      <- file.path(project_root, "outputs")

# Load dataset
data_long <- read_excel(file.path(intermediate_dir, "eqosogi_score.xlsx"))

# Build feature set robustly
features_df <- data_long %>%
  group_by(country) %>%
  arrange(year) %>%
  nest() %>%
  mutate(
    slope_total = map_dbl(data, ~ coef(lm(es_ ~ year, data = .x))[2]),  # Total slope
    delta = map_dbl(data, ~ last(.x$es_) - first(.x$es_)),              # Total change
    first_value = map_dbl(data, ~ first(.x$es_)),                       # Initial value
    last_value = map_dbl(data, ~ last(.x$es_)),                         # Final value
    sd = map_dbl(data, ~ sd(.x$es_)),                                   # Standard deviation
    mean = map_dbl(data, ~ mean(.x$es_)),                               # Mean
    max_jump = map_dbl(data, ~ max(diff(.x$es_))),                      # Largest increase
    year_of_max_jump = map_dbl(data, ~ .x$year[which.max(diff(.x$es_)) + 1]),  # Year of largest jump
    min_jump = map_dbl(data, ~ min(diff(.x$es_))),                      # Largest decrease
    n_neg_jumps = map_dbl(data, ~ sum(diff(.x$es_) < 0)),               # Number of declines
    prop_zero_jumps = map_dbl(data, ~ mean(diff(.x$es_) == 0)),         # % with no change
    slope_early = map_dbl(data, ~ coef(lm(es_ ~ year, data = filter(.x, year >= 1960 & year < 1981)))[2]),
    slope_mid   = map_dbl(data, ~ coef(lm(es_ ~ year, data = filter(.x, year >= 1981 & year < 2001)))[2]),
    slope_late  = map_dbl(data, ~ coef(lm(es_ ~ year, data = filter(.x, year >= 2001)))[2]),
    delta_mid   = map_dbl(data, ~ .x$es_[.x$year == 2000] - .x$es_[.x$year == 1981]),
    delta_late  = map_dbl(data, ~ .x$es_[.x$year == 2024] - .x$es_[.x$year == 2001]),
    iqr = map_dbl(data, ~ IQR(.x$es_))
  ) %>%
  select(-data) %>%
  ungroup()

# ---- function: from a feature vector -> a one-row data.frame with metrics ----
cluster_metrics_df <- function(features_pick, k = 4, seed = 123) {
  f <- features_df %>%
    dplyr::select(country, dplyr::any_of(features_pick)) %>%
    tidyr::drop_na()
  if (nrow(f) < k) stop("Not enough countries with complete data")
  
  X <- scale(dplyr::select(f, -country))
  set.seed(seed)
  km <- kmeans(X, centers = k)
  
  mean_sil <- NA_real_
  if (length(unique(km$cluster)) >= 2) {
    sil <- cluster::silhouette(km$cluster, dist(X))
    mean_sil <- mean(sil[, 3])
  }
  
  db <- clusterCrit::intCriteria(as.matrix(X), as.integer(km$cluster), "Davies_Bouldin")$davies_bouldin
  ch <- clusterCrit::intCriteria(as.matrix(X), as.integer(km$cluster), "Calinski_Harabasz")$calinski_harabasz
  
  data.frame(
    "Mean Silhouette Score"   = sprintf("%.3f", mean_sil),
    "Davies-Bouldin Index"    = sprintf("%.3f", db),
    "Calinski-Harabasz Index" = sprintf("%.3f", ch),
    check.names = FALSE
  )
}

# ------------------------------------------------------------------------------- #
# Table A2: Internal validation metrics for alternative clustering specifications #
# ------------------------------------------------------------------------------- #

## Number of clusters: 4

# Baseline: All features
features_pick <- c("slope_total","delta","first_value","last_value","sd","mean","max_jump",
                   "year_of_max_jump","min_jump","n_neg_jumps","prop_zero_jumps","slope_early",
                   "slope_mid","slope_late","delta_mid","delta_late","iqr")
metrics_all <- cluster_metrics_df(features_pick)  # one row with labels
features_pick <- c("slope_total","delta","first_value","last_value","sd","mean","max_jump")
metrics_set1 <- cluster_metrics_df(features_pick)  # one row with labels
features_pick <- c("slope_total","delta","last_value","sd","mean","max_jump")
metrics_set2 <- cluster_metrics_df(features_pick)  # one row with labels
features_pick <- c("slope_total","delta","last_value","sd","mean","max_jump","delta_late")
metrics_set3 <- cluster_metrics_df(features_pick)  # one row with labels
features_pick <- c("slope_total","delta","last_value","sd","mean","max_jump","slope_late","delta_late")
metrics_set4 <- cluster_metrics_df(features_pick)  # one row with labels

table_k4 <- rbind(metrics_all, metrics_set1, metrics_set2, metrics_set3, metrics_set4)

table_k3 <- read_excel(file.path(intermediate_dir, "TableA2_k3.xlsx"))

# table_k3 y table_k4 ya existen y tienen 5 filas y 3 columnas:
# "Mean Silhouette Score", "Davies-Bouldin Index", "Calinski-Harabasz Index"

# Etiquetas de fila (columna "Features"), en el mismo orden de las filas 1–5
features <- c(
  "Baseline: All features",
  "Set 1: slope_total, delta, first_value, last_value, sd, mean, max_jump",
  "Set 2: slope_total, delta, last_value, sd, mean, max_jump",
  "Set 3: slope_total, delta, last_value, sd, mean, max_jump, delta_late",
  "Set 4: slope_total, delta, last_value, sd, mean, max_jump, slope_late, delta_late"
)

# (opcional) chequeo rápido
stopifnot(length(features) == nrow(table_k3),
          nrow(table_k3) == nrow(table_k4))

# Nombres de las métricas
metric_cols <- names(table_k3)

# Construir columnas intercaladas k = 3 / k = 4
body_A2 <- map_dfc(metric_cols, function(col) {
  tibble(
    !!paste0(col, " (k = 3)") := table_k3[[col]],
    !!paste0(col, " (k = 4)") := table_k4[[col]]
  )
})

# Armar tabla final con columna "Features"
tableA2 <- tibble(Features = features) %>%
  bind_cols(body_A2)

# Export to Excel
write_xlsx(tableA2, paste0(outputs_dir, "/TableA2.xlsx"))
